AFOSR-tfc  9  1  06  9.1 


AD-A238  158 

IIIIIHIHII 

A  STUDY  OF  THE  BEHAVIOR  AND  MICROMECHANICAL  MODELLING  OF 

GRANULAR  SOIL 


DTIC 


ELECTE 
JUL1  2  1991 


A*. 


VOLUME  in 

A  NUMERICAL  INVESTIGATION 
OF  THE  BEHAVIOR  OF  GRANULAR  MEDIA 
USING  NONLINEAR 
DISCRETE  ELEMENT  SIMULATIONS 

by 

Emmanuel  Petrakis 
Ricardo  Dobry 
Tang-Tat  Ng 

.release*  Li  Liu 

Kited. 


Prepared  under  Contract  No.  AFOSR-89-0350 
United  States  Air  Force 
Office  of  Scientific  Research 
Bolling  Air  Force  Base 


r—  ;._j 

fcj  o 

h- , ,  * 
6  £’ 


0) 

* _ »  , 


Department  of  Civil  Engineering 
Rensselaer  Polytechnic  Institute 

Troy,  NY  12180-3590  ,v  ... 

*• 

J-*5 

May  1991 


-DfSfRlWTIOH  VT.'VTggLi 

Approved  for  public  release: 
Distribution  Unlimited^ — 


91-04782 

IH1IRHUHH 


91  7  H 


076 


REPORT  DOCUMENTATION  PAGE 

FO"*  A00rov«0 

OMB  No.  0704-014B 

f  '■  n,!>  7i. i1' 'L'.'::1. 

May  22,  1991  Final  1/6/89-5/15/91 

4b  mu  AMO  SURTTTU 

A  Study  of  the  Behavior  and  Micromechanical  Modelling  of 
Granular  Soil  , f  i  ^ 

VolovYiP  3. 

5.  PUNOING  NUMRIRS 

Grant  AF0SR-89-0350 

PR  2302/Cl 

1  AUTHORS) 

Emmanuel  Petrakis,  Ricardo  Dobry,  Paul  Van  Laak,  Panos 
Kotsanopoulos,  Tang-Tat  Ng,  and  Li  Liu 

7.  PCRPORMMG  ORGANIZATION  NAMI(S)  ANO  AOORISSilS) 

Civil  Engineering  Department 

Rensselaer  Polytechnic  Institute 

Troy,  NY  12180-3590 

«.  PIRPORMtNG  ORGANIZATION 

RIPORT  NUMRIR 

t.  SPONSORING/  MONITORING  AGINCY  NAM«S)  ANO  AOORISKIS) 

AFOSR/NA 

Bldg.  410 

Bolling  AFB 

Washington,  DC  20332-6448 

10.  SPONSORING/  MONITORING 

AGINCY  RIPOKT  NUMRIR 

ffpOZK'- 

SUPPllMIMTAKY  NOUS 

tk  OtSiRMNiTlOM /  AVARJkRRjTY  STAUDMNT 

Approved  for  Public  Release;  Distribution  Unlimited 

13.  ARSTRACT  (Mtjomvm  2QQ  worGN 

A  comprehensive  research  effort  has  been  conducted  on  constitutive 
and  micromechanical  modelling  of  granular  soil.  This  Includes:  1)  the  development 
of  a  new  constitutive  relation  for  granular  media  based  on  the  contact  law  between 
two  spheres;  2)  an  experimental  Investigation  on  the  stress-strain  response  of 
a  glass  bead  matetlal  with  46  monotonic  and  cyclic  experiments  on  hollow  cylinder 
specimens,  most  o'  them  constant  mean  stress  tests  to  measure  devlatoric  response 
and  behavior  of  iltlal  and  subsequent  yield  loci;  and  3)  numerical  simulations 
of  the  behavior  r  '  granular  media  using  the  discrete  element  method. 

The  proposed  constitutive  law  captures  a  number  of  key  aspects  of 
the  observed  stress-strain  behavior  of  granular  soils,  and  It  predicts  well  the 
experiments  on  glass  beads.  Novel  aspects  of  the  proposed  model  Include  yield 
cones  parallel  to  the  failure  envelope,  and  a  basic  relation  between  the  field 
of  elastoplastlc  moduli  and  the  elastic  constants  of  the  material. 
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Micromechanics;  Constitutive  Law;  Granular  Media;  Sand;  Yield 
Surface  Distortion;  Hollow  Cylinder  Experiment;  Discrete  Element 
Simulations;  Plasticity 
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The  report  consists  of  three  volumes,  as  follows: 

Volume  I:  "A  Constitutive  Law  for  Granular  Materials  Based  on  the  Contact  Law 
Between  Two  Spheres,"  by  R.  Dobry  and  E.  Petra kis. 

Volume  II:  "An  Experimental  Investigation  of  the  Behavior  of  Granular  Media  Under 
Load,"  by  E.  Petrakis,  R.  Dobry,  P.  Van  Laak,  and  P.  Kotsanopoulos. 

Volume  III:  “A  Numerical  Investigation  of  the  Behavior  of  Granular  Media  Using 
Nonlinear  Discrete  Element  Simulations,"  by  E.  Petrakis,  R.  Dobry,  T.-T.  Ng,  and 
L.  Liu. 
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ABSTRACT 

This  report  presents  the  results  of  a  study  on  the  monotonic  and  cyclic  stress-strain  behavior  of 
quartz  sands  and  other  granular  media.  Random  arrays  of  elastic  rough  spheres  are  used  to 
represent  the  medium  in  the  numerical  simulations  which  utilize  one,  two  and  three  sizes  of 
spheres.  In  all  cases  the  spherical  grains  are  assigned  the  properties  of  quartz.  Results  from  a 
2D  array  of  531  grains  and  from  a  3D  array  of  395  particles  are  displayed  and  discussed. 

A  rigorous  representation  of  granular  stress-strain  response  at  small  strains  and  during  cyclic 
loading  requires  the  correct  modelling  of  the  force-displacement  behavior  at  the  contact  between 
two  grains.  This  was  achieved  through  a  numerical  solution  of  the  (Hertz-Mindlin)  contact 
problem  between  two  spheres  developed  by  the  authors  which  uses  the  theory  of  plasticity  and 
kinematic  hardening  to  describe  the  phenomenon.  The  discrete  element  finite  difference  method 
proposed  by  Cundall  and  Strack  is  used  to  calculate  the  interactions  between  particles  including 
geometric  rearrangement  and  formation  of  new  contacts.  A  periodic  space  solution  is 
implemented  to  avoid  undesirable  boundary  effects.  All  these  aspects  were  incorporated  by  the 
authors  into  program  CONBAL-P  (CONtact  truBAL-Parallel).  Because  of  the  large  number  of 
contacts,  the  complex  nonlinear  character  of  the  contact  law,  and  the  explicit  nature  of  the  finite 
difference  scheme,  a  supercomputer  is  required  to  run  CONBAL-P,  which  has  ►  een  optimized 
for  a  vector  and  parallel  environment. 

The  results  discussed  in  the  paper  focus  on  identifying  the  micromechanical  phenomena 
responsible  for  the  motion  and  distortion  of  the  yield  loci  of  a  granular  medium  in  stress  space 
during  shear  loading  that  were  observed  in  the  laboratory  experiments  of  Volume  II.  The  array 
is  first  isotropically  compressed  and  then  is  sheared  at  constant  mean  stress.  The  shear  loading 
consists  of  monotonically  prestraining  the  medium  to  a  predetermined  level  which  is  followed 
by  unloading  and  probing  in  different  directions  to  establish  the  current  size  and  shape  of  a  yield 
locus.  It  has  been  found  that  the  yield  locus  distorts  in  the  direction  of  loading  in  a  way  similar 
to  that  observed  in  the  laboratory  experiments  on  glass  beads.  This  is  in  contrast  to  the  current 
assumption  for  soils  that  the  yield  locus  translates  without  distortion.  These  results  are  relevant 
to  the  formulation  of  constitutive  relation  for  granular  soil  since  they  indicate  that  the  anisotropic 
fabric  which  develops  during  loading  needs  to  be  modelled  in  order  to  capture  the  macroscopic 
phenomena.  The  simulations  also  provide  a  wealth  of  statistical  information  which  is  used  to 
correlate  the  macroscopic  stress-strain  response  observed  in  the  laboratory  with  micromechanical 
phenomena  at  the  particle  contact  level. 
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1.  INTRODUCTION 
1.1  General 

Over  the  past  30  years  considerable  attention  has  been  given  to  the  development  of  constitutive 
laws  for  engineering  materials  (Hill  1950,  Prager  1955,  Mroz  1967,  Dafalias  and  Popov  1976, 
Drucker  and  Palgen  1981,  Yen  and  Eisenberg  1987).  Among  other  formulations,  the  existing 
models  are  based  on  the  theories  of  elasticity,  hypoelasticity,  plasticity  and  viscoplasticity. 
Despite  the  large  number  of  models,  there  is  no  consensus  yet  within  the  scientific  research 
community  on  the  best  approach.  However,  the  models  based  on  the  theory  of  plasticity  or 
viscoplasticity  appear  to  be  most  promising.  Especially  critical  is  the  need  for  a  constitutive 
relation  for  soils  required  for  analyses  of  onshore  and  offshore  structures  and  their  foundations, 
as  well  as  earth  structures,  subjected  to  a  variety  of  natural  and  man-made  static  and  dynamic 
loading  environments.  Soil  behavior  is  more  complex  than  that  of  other  materials  due  to  its 
particulate  nature,  which  makes  it  pressure-dependent  and  nonlinear  inelastic  even  at  small 
strains.  A  large  number  of  models  has  been  proposed  for  soils  based  on  the  theory  of  plasticity, 
including  those  of  DiMaggio  and  Sandler  (1971),  Baladi  and  Rohani  (1979),  Lade  (1977), 
Prevost  (1978,  1985),  Dafalias  and  Herrmann  (1982).  In  these  models,  the  total  strain  is  equal  to 
the  sum  of  the  elastic  and  plastic  strain  increments,  de^  =  dr,,-  +  de*^  with  these  increments  being 
rate  independent.  A  variety  of  associative  and  non-associative  flow  rules  have  been  proposed  for 
the  plastic  strain  increment,  of  die  form: 

<u5.<u-i£-  a) 

dOy 

where  A.  is  a  coefficient  of  proportionality  and  g(oij)  is  the  plastic  potential  function,  which  may 
or  may  not  coincide  with  the  yield  function,  f(oy),  at  which  plastic  strains  develop.  It  is  typically 
assumed  that  these  yield  surfaces  (loci)  translate  in  stress  space  without  distortion  (kinematic 
strain  hardening),  sometimes  increasing  or  decreasing  in  size  (isotropic  strain  hardening). 

The  limitation  of  these  models  is  that  they  are  phenomenological.  They  have  been  typically 
developed  from  a  manageable  mathematical  formulation,  and  they  have  been  calibrated  and 
modified  by  interpreting  macroscopic  experimental  results,  while  ignoring  the  underlying 
micromechanical  phenomena.  As  a  result,  the  existing  plasticity  models  for  soils  are  in  need  for 
constant  refinement  when  needed  for  cases  very  different  from  the  one  the  model  was  originally 
developed  and  calibrated  for. 
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The  current  situation  in  metal  plasticity  is  quite  different  Although  the  modelling  of  the 
nonlinear  behavior  of  metals  started  on  a  similar  phenomenological  basis,  there  has  been  a  shift 
in  the  last  20  years  or  so  toward  formulating  the  metal  response  with  due  consideration  of 
micromechanical  principles  (Budiansky  and  Wu  1962,  Lin  and  Ito  1965,  1966).  Recently,  this 
has  been  enhanced  by  specific  experiments  and  micromechanical  (transmission  electron 
microscopy  -  TEM)  measurements  (Stout  et  al.  1985,  Helling  et  al.  1986).  The  situation  is 
analogous  in  the  modelling  of  more  complex  composite  materials,  where  experiments  and 
micromechanical  analytical  simulations  are  combined  to  create  the  corresponding  constitutive 
law  (Dvorak  1987,  Dvorak  et  al.  1988). 

A  similar  effort  in  soils  has  begun  at  RPI  by  the  authors,  in  which  a  new  constitutive  law  for 
granular  media  is  being  developed  which  evolves  from  several  micromechanical  studies, 
including  numerical  simulations  of  2-D  and  3-D  random  arrays  of  spheres  and  laboratory 
experiments  on  glass  beads  (Petrakis  and  Dobry  1989).  This  paper  presents  the  results  of  2  and 
3-D  simulations  of  the  mechanical  behavior  of  random  assemblages  of  rough,  elastic  spheres. 

1.2  Studies  on  Polvcrvstalline  Aggregates 

Starting  in  the  1950's,  researchers  have  simulated  the  elastic-plastic,  stress-strain  behavior  of 
polycrystalline  aggregates  through  analytical,  semi-analytical,  and  numerical  micromechanical 
techniques  (Hershey  1954,  Budiansky  and  Wu  1962,  Lin  and  Ito  1965,  1966,  Hill  1967,  Canova 
et  al.  1985).  This  micromechanical  work  has  not  supported  the  continuum  mechanics 
hypotheses  of  either  pure  kinematic  or  isotropic  hardening  behavior,  but  has  predicted  a 
combined  translation  (kinematic  hardening)  and  distortion  of  the  yield  surfaces  in  the  direction 
of  loading. 

The  micromechanical  approach  commonly  used  to  analyze  the  elastic-plastic  behavior  of 
polycrystals  assumes  that  they  are  an  assemblage  of  equal  anisotropic  monocrystals  (single 
crystals)  randomly  oriented  in  space.  This  results  in  an  isotropic  polycrystal  if  the  spatial 
distribution  of  the  orientations  is  statistically  uniform.  A  monocrystal  has  n  sliding  planes  with 
each  plane  having  m  sliding  directions,  and  with  every  sliding  direction  corresponding  to  a  pair 
of  parallel  yield  planes  in  stress  space.  In  the  limiting  case  in  which  an  infinite  number  of 
possible  crystal  orientations  is  assumed,  this  infinitely  sided  polyhedron  becomes  a  curved  yield 
surface. 
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Plastic  strain  in  the  aggregate  is  caused  by  sliding  of  one  of  the  sliding  planes  occurring  in  a 
family  of  similarly  oriented  crystals.  After  sliding  has  occurred  in  a  number  of  these  families, 
each  surface  of  the  polyhedron  mentioned  above  expands  and  shifts  differently.  These  sliding 
directions  are  all  more  or  less  parallel  to  the  direction  of  the  plane  of  the  maximum  shear  stress 
acting  on  the  aggregate.  As  the  aggregate  is  loaded  further  beyond  the  elastic  range,  more 
crystals  and  crystal  families  slide,  and  increasingly  more  yield  planes  pass  through  the  loading 
point  on  the  yield  surface.  The  yield  planes  of  different  orientations  intersect  at  that  point  on  the 
yield  surface  and  form  a  comer  or  vertex. 

This  vertex  is  not  easily  observed  during  testing.  One  reason  is  that  the  very  large  number  of 
monocrystal  orientations  smooths  the  effect,  which  appears  as  a  "smooth  vertex"  or  distortion  of 
the  yield  surface  in  the  direction  of  loading,  rather  than  a  sharp  comer.  This  distortion  of  the 
yield  surface  associated  with  the  vertex  reflects  the  "memory"  the  material  has  of  prestraining  in 
the  direction  of  loading.  The  existence  of  this  yield  surface  distortion  in  metals  has  been 
observed  experimentally  by  a  number  of  researchers  in  several  polycrystalline  aggregates, 
including  aluminum  alloys,  brass  and  magnesium  (Naghdi  et  al.  1958,  Phillips  1968,  Kelley  and 
Hosford  1968,  Phillips  et  al.  1970,  Shiratori  et  al.  1976,  Helling  et  al.  1986). 

The  late  Professor  Phillips  developed  a  testing  procedure  to  seek  and  compute  the  initial  and 
subsequent  yield  surfaces  of  aluminum  and  their  motion  in  stress  space  (Fig.  la).  This 
procedure  is  now  widely  used  in  experimental  plasticity  studies  of  metals  and  metal  matrix 
composites  (Rousset  1985,  Stout  et  al.  1985,  Dvorak  1987,  Dvorak  et  al.  1988).  The  experiment 
are  typically  conducted  by  applying  a  combination  of  tensile  and  torsional  shear  stresses  to  a 
hollow  cylindrical  specimen.  In  these  tests,  the  loading  stops  and  reverses  as  soon  as  a  point  on 
the  yield  surface  is  reached,  as  defined  by  a  certain  deviation  from  the  linear  portion  of  the 
stress-strain  curve.  Then  the  small  strain  loading  (probe)  continues  until  another  point  on  the 
surface  is  reached,  the  load  is  then  reversed,  and  then  increased  until  enough  points  have  been 
obtained  to  accurately  define  the  yield  surface  (locus).  Fig.  lb  clearly  shows  the  characteristic 
distortion  of  the  yield  surface  in  aluminum  in  the  direction  of  loading.  While  the  yield  surface 
(for  a  given  temperature  )  in  the  t-<J,  space  is  an  ellipse,  the  subsequent  yield  surfaces  have 
distorted  and  become  pointed  in  the  direction  of  loading  (a),  while  becoming  flatter  in  the 
opposite  direction  (b).  As  a  result,  the  size  of  the  yield  surface  shrinks  in  the  direction  of 
loading  while  staying  constant  in  the  other  direction. 

Laboratory  results  such  as  these  have  made  possible  the  Unking  of  the  micromechanical  theory 
with  experiments  (Stout  et  al.  1985,  Helling  et  al.  1986),  and  have  led  to  a  new  family  of 
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constitutive  laws  (Phillips  and  Weng  1975,  Eisenberg  and  Yen  1981,  1984,  Yen  and  Eisenberg 
1987)  which  incorporate  the  above  findings. 

1.3  Modelling  of  Granular  Soil 

The  behavior  of  a  sand  aggregate  is  similar  to  that  of  a  polycrystal,  since  the  individual  groups 
of  grains  or  grain  packings  within  the  sand  may  be  considered  in  first  approximation  to  behave 
like  randomly  oriented  crystals.  The  main  difference  is  that  the  properties  of  these  packings  are 
now  pressure  dependent.  For  example,  a  simple  cubic  array  of  equal  spheres  is  a  pressure 
dependent  monocrystal  having  three  sliding  planes  (n=3),  with  each  sliding  plane  containing  two 
sliding  directions  90°  apart  (m=2). 

As  in  the  case  of  polycrystalline  aggregates,  each  sliding  plane  in  each  of  the  packings 
corresponds  to  a  pair  of  parallel  yield  planes  in  stress  space.  The  macroscopic  yield  surface  of 
the  material  is  the  surface  bounding  the  yield  polyhedron,  the  sides  of  which  are  formed  by  the 
intersection  of  these  yield  planes.  Plastic  strain,  defined  as  irrecoverable  deformation,  is  the 
result  of  a  slide  in  at  least  one  of  these  packings.  This  is  directly  analogous  to  the  definition  of 
yielding  in  metals  where  plastic  strain  is  caused  by  a  slip  in  one  of  the  slip  systems  and  the  term 
"sliding"  in  soils  is  the  equivalent  of  "slip"  in  metals.  However,  soils  exhibit  nonlinear  inelastic 
stress-strain  behavior  even  at  small  strains  which  is  the  result  of  the  force-deformation 
nonlinearity  at  the  interparticle  contacts.  Nevertheless,  granular  soils  experience  little  or  no 
sliding  between  particle  contacts  and  thus  exhibit  nondestructive  behavior  and  no  dilation  up  to 
the  so  called  "threshold  shear  strain",  yt,  (1X104  <yt<  2X104 ;  Dobry  et  al.  1982).  Therefore,  in 
granular  soils,  loading  below  the  threshold  has  some  important  features  in  common  with  elastic 
loading  in  metals  before  the  initial  yield  surface  is  reached.  In  soils  the  loading  below  y,  though 
nondestructive,  is  inelastic  and  does  include  some  plastic  yielding  due  to  localized  slipping 
within  the  intergranular  contact  areas.  If  this  localized  slipping  effect  below  yt  is  neglected  in 
first  approximation,  one  possible  definition  of  an  initial  yield  surface  in  cohesionless  soils  could 
be  the  surface  in  stress  space  where  y,  is  reached,  that  is  approximately  at  strains  of  less  than 
2X1(H 
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2.  NUMERICAL  METHODS  AND  PROGRAM  CONBAL-P 
2.1  General 

The  development  of  a  realistic  computer  program  to  simulate  granular  media  by  random  arrays 
of  elastic,  rough  particles  presupposes  a  number  of  steps  and  decisions  related  to:  shapes,  sizes, 
properties  and  initial  spatial  distributions  of  grains,  boundary  conditions,  force-deformation  law 
between  two  adjacent  particles  in  contact  A  numerical  scheme  needs  to  be  chosen  which  will 
allow  for  individual  grain  interactions  and  rearrangements  including  disappearance  of  old 
contacts  and  formation  of  new  ones.  In  what  follows,  some  of  the  solutions  provided  by  the 
authors  in  development  of  their  method  and  associated  computer  code  are  discussed. 

The  problem  of  the  contact  of  two  elastic,  elliptical  semi-infinite  bodies  subjected  to  a  normal 
force  was  first  studied  by  Hertz  (1882).  Hertz  demonstrated  that  the  normal  force-deformation 
behavior  at  the  contact  is  nonlinear  elastic.  Cattaneo  (1938),  Mindlin  (1949)  and  Mindlin  and 
Deresiewicz  (1953)  addressed  the  problem  of  the  contact  of  two  identical  elastic,  rough  spheres 
subjected  to  a  combination  of  normal  and  tangential  forces  and  presented  a  number  of  closed 
form  solutions  for  each  of  several  specific  load  histories. 

Seridi  and  Dobry  (1984)  and  Dobry  et  al.  (1991)  developed  a  numerical  solution  for  the 
force-deformation  behavior  of  two  elastic,  rough  spheres  in  contact  under  a  combination  of 
arbitrarily  varying  normal  and  tangential  forces,  which  was  then  coded  as  program  CONTACT. 
This  model  is  based  on  the  incremental  theory  of  plasticity,  uses  an  infinite  number  of  yield 
surfaces  (which  correspond  to  an  infinite  number  of  contact  annuli)  and  assumes  kinematic 
hardening  (Fig.  2).  Therefore,  its  main  features  are  very  similar  to  those  of  the  plasticity 
stress-strain  models  for  engineering  materials  described  previously. 

The  numerical  solution  described  previously,  implemented  in  computer  program  CONTACT 
was  subsequently  tested  to  verify  that  it  reproduces  accurately  the  analytical  solutions  obtained 
by  Mindlin  and  Deresiewicz  (1953).  Once  this  was  done,  program  CONTACT  was  made  a 
subroutine  which  could  be  used  in  Finite  Element  and  Discrete  Element  (Cundall  and  Strack 
1979)  simulations. 
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2.2  The  Discrete  Element  Method 

Many  studies  have  been  performed  using  particulate  models  to  interpret  and  model  the  behavior 
of  cohesionless  soils  and  other  granular  materials.  Most  of  these  investigations  have  been  also 
included  measurements  in  actual  granular  soils,  as  well  as  in  regular  or  random  arrays  of  spheres 
(3-D)  or  disk/rod  assemblies  (2-D);  a  number  of  them  have  dealt  with  the  load-deformation 
characteristics  at  the  contact  between  two  elastic  bodies  possessing  friction. 

Of  special  interest  are  the  investigations  which  have  studied  the  stress-strain  behavior  of  granular 
arrays  considering  the  elasticity  of  the  particles  and  the  corresponding  compliances  at  the 
contacts  (Duffy  and  Mindlin  1958,  Duffy  1959,  Deresiewicz  1958,  Walton  1987,  Petrakis  and 
Dobry).  Serrano  and  Rodrigues-Ortiz  (1973)  suggested  a  method  for  generating  a  random 
configuration  of  unequal  disks  or  spheres  having  a  prescribed  grain  size  distribution.  Cundall 
and  Strack  (1979)  used  a  similar  approach  in  conjunction  with  their  distinct  element  method  to 
successfully  simulate  the  mechanical  behavior  of  arrays  of  disks  and  spheres  under  a  variety  of 
loading  conditions.  In  their  method,  an  explicit  finite  difference  formulation  is  used  to 
determine  the  static  response  of  the  array  to  applied  strains  (program  TRUBAL)  or  boundary 
displacement  (program  BALL). 

Program  TRUBAL  uses  a  periodic  space  in  a  form  of  a  cube  (3-D)  or  square  (2-D),  to  minimize 
the  effect  of  boundaries  and  allow  a  relatively  small  number  of  particles.  A  uniform  strain  field 
is  applied  to  all  particles  and  the  corresponding  contact  forces  are  calculated  from  the  relative 
displacements  between  neighboring  spheres.  As  illustrated  by  Fig.  3,  these  contact  forces,  Pit 
produce  in  each  spheres  a  resultant  unbalanced  force  ZPit  and  unbalanced  moment  ZM;,  which 
induce  linear  and  angular  accelerations  to  the  particle.  These  accelerations  are  used  in  turn  to 
compute  new  particle  position  at  the  end  of  a  time  increment,  t,  and  new  contact  forces.  The 
process  is  repeated  until  static  equilibrium  is  achieved  (  ZP;  =  ZMS  =  0).  The  method  has  the 
advantage  of  decreasing  substantially  the  required  computer  memory,  as  no  large  stiffness 
matrices  need  to  be  calculated  and  inverted.  However,  the  execution  time  is  large  due  to  the 
great  number  of  iterations  needed  to  assure  static  equilibrium.  Initially,  TRUBAL  was  using  a 
linear,  non-pressure  dependent  law  to  describe  the  force-deformation  behavior  at  the  interparticle 
contacts;  a  later  version  of  the  program  (Zhang  and  Cundall,  1986)  used  a  linear, 
pressure-dependent  force-deformation  relationship. 

The  existing  force-deformation  law  at  the  interparticle  contacts  in  program  TRUBAL  was 
modified  by  the  rigorous  numerical  solution  to  the  contact  problem  (subroutine  CONTACT).  To 
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do  this,  a  new  array  was  introduced  in  TRUBAL  which  stores  the  load  history  at  each  contact 
(information  on  the  yield  surfaces,  elastoplastic  moduli,  position  of  apexes,  force  point  etc).  The 
addition  of  this  complex  and  computationally  intensive  nonlinear  force-deformation  law  at  each 
interparticle  contact  increased  the  computational  demands  of  the  program  which  could  no  longer 
be  feasible  to  run  on  a  regular  computer.  A  preliminary  study  shows  that  by  adding  this 
complex  force-deformation  law  to  the  program,  the  computation  time  increases  by  a  factor  of 
six. 

As  a  result  of  this,  the  program  was  further  modified  by  vectorizing  the  code  to  run  efficiently 
on  a  supercomputer  with  vector  facilities.  Recently  the  code  has  been  parallelized  as  well,  so  as 
to  take  advantage  of  a  computing  environment  of  six  processors.  As  shown  in  Fig.  3,  since  each 
sphere  is  treated  individually,  aifferent  contact  forces  between  two  different  spheres  can  be 
calculated  simultaneously.  To  achieve  this,  the  specimen  is  divided  into  different  zones  and  the 
contact  forces  between  different  spheres  in  different  zones  are  calculated  at  the  same  times. 
Calculations  in  each  zone  is  assigned  as  a  sub  task  for  each  of  the  six  processors  of  the 
supercomputer  used  in  this  study  (IBM  3090-600E).  For  a  loading  simulation  of  a  very  small 
3-D  specimen  (104  spheres),  the  parallel  code  in  the  six  processor  environment  of  the  Cornell 
National  Supercomputer  Facility  runs  more  than  twice  as  fast  (in  real  clock  time)  as  it  would  run 
if  only  one  processor  was  used.  The  CPU  time  increases,  as  expected  (see  Table  1). 


TABLE  1 


Processing  Mode 

CPU  time 

Clock  time 

Serial  (vector) 

141  min 

210  min 

Parallel  (6  processors) 

140  min 

60  min 

This  performance  increase  is  more  dramatic  and  the  computational  efficiency  improves  when 
larger  specimens  (more  spheres)  are  used.  This  new  modified  program  is  named  CONBAL-P 
which  is  essentially  a  version  of  TRUBAL  incorporating  CONTACT  and  optimized  for  the 
supercomputer.  Therefore,  CONBAL  is  very  similar  to  original  program  TRUBAL,  except  that 
now  the  effect  of  the  normal  force  at  the  intergranular  contacts  is  included  and  rigorously 
modelled.  The  program  CONBAL-P  can  run  in  2-D,  in  which  case  the  behavior  of  a  monolayer 
of  spheres  whose  centers  lie  on  the  same  plane  is  analyzed,  or  in  3-D  where  the  behavior  of  3-D 
assemblies  of  spheres  is  modelled.  Program  CONBAL-P  has  been  requiring  up  to  150MB  of 
primary  memory  storage  to  run  (depending  on  the  complexity  of  the  loading  case  and  the 
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complexity  of  the  law  modelling  the  contact  behavior)  and  the  computation  time  for  a  typical  job 
is  6  CPU  hours  (vector)  on  an  IBM-3090-600E  for  a  2-D  simulation  and  18  CPU  hours  for  a  3-D 
simulation  utilizing  the  full  Mindlin  solution  without  considering  particle  rotation. 

Since  the  interparticle  force-deformation  relationships  developed  by  Cattaneo  (1938),  Mindlin 
(1949)  and  Mindlin  and  Deresiewicz  (1953)  apply  only  to  spheres  of  the  same  size  and 
properties,  only  random  arrays  of  equal  and  moderately  unequal  spheres  were  developed  and 
studied.  Subroutine  CONTACT,  and  thus  CONBAL-P  can  in  principle  analyze  only  random 
arrays  of  equal  spheres;  however,  the  same  program,  with  certain  assumptions,  can  handle  in 
first  approximation  moderately  unequal  spheres  of  the  same  material.  Moreover,  the  same 
analytical  solutions  do  not  consider  rolling  of  the  spheres;  consequently  rolling  is  not  permitted 
in  the  present  version  of  CONBAL-P  where  the  fully  nonlinear  contact  model  is  utilized.  In  the 
case  that  rotation  is  not  allowed,  the  coefficient  of  friction  of  the  material  of  the  spheres  has  been 
reduced  accordingly  (Rowe  1962).  A  second  version  of  the  program  allows  for  particle  rotation 
and  a  simplified  contact  model  is  used,  which  assumes  only  one  diameter  of  the  contact  annulus 
(sliding  -  no  sliding  condition  using  one  yield  surface  in  subroutine  CONTACT).  This  program 
takes  considerably  less  time  to  run  and  has  been  running  on  workstations  as  well. 

The  accuracy  of  CONBAL-P  has  been  checked  by  comparing  its  results  at  both  small  and  large 
strains  to  a  number  of  available  analytical  solutions  with  excellent  agreement  for  both  static  and 
cyclic  loading  under  a  number  of  boundary  conditions. 

The  simulations  using  CONBAL-P  can  be  of  great  importance  to  granular  media  research, 
because  the  provide  detailed  information  on  any  micromechanical  variables,  such  as  the 
distribution  of  contact  angles,  the  distribution  of  the  mobilized  angle  (the  angle  between  the 
shear  force  and  contact  normal),  the  average  number  of  contacts  per  sphere  (coordination 
number,  CN)  and  their  evolution  with  loading.  These  simulations  can  give  the  same  information 
that  Magnetic  Resonance  Imaging  (MRI)  would  give  if  applied  to  a  laboratory  specimen.  The 
micromechanical  information  could  then  be  used  to  interpret  the  macroscopic  behavior. 


11 


3.  NUMERICAL  SIMULATIONS 
3.1  Monotonic  Loading  Simulations 

As  a  first  step,  it  was  decided  to  start  simulating  phenomena  in  the  small  strain  range  (10-6  <yt< 
10-3),  and  at  a  later  stage  focus  on  the  large  strain  (yt  >  10  3),  fully  nonlinear  inelastic 
macroscopic  behavior.  Moreover,  it  was  also  decided  to  compare  the  results  of  CONBAL  with 
answers  obtained  using  analytical  (Self  Consistent,  Petrakis  and  Dobry  1991)  and  numerical 
(nonlinear  Finite  Element,  Petrakis  and  Dobry  1991)  procedures  developed,  as  well  as  with 
experimental  data  in  the  literature  and  on  glass  beads  presented  in  Volume  n.  Results  of  small 
strain  parametric  studies  on  two  2-D  random  arrays  of  spheres,  one  with  477  equal  particles  and 
area  porosity  (calculated  on  the  plane  passing  through  the  sphere  centers)  n  =  0.154  and  the  other 
with  531  particles  with  a  ratio  of  radii  Rj/  Rj  =  1.5  and  area  porosity  n  =  0.182  appear  in  Petrakis 
et  al.  (1989).  Parametric  studies  on  a  3-D  array  of  104  particles  with  R,/  R^  =  1.5  and  porosity  n 
=  0.359  appear  in  Ng  and  Dobry  (1990). 

The  geometrical  configuration  of  the  477  equal  particle  array  subjected  to  an  isotropic  confining 
pressure  o0=  91  KPa  is  shown  in  Fig.  4.  The  configuration  of  the  531-particle  array  under  o0  = 
132  KPa  is  shown  in  Fig.  5.  In  both  figures  the  circles  represent  the  spheres  and  the  rectangles 
the  relative  magnitude  and  direction  of  the  contact  force.  There  are  four  different  rectangle 
widths,  with  each  one  of  them  corresponding  to  a  range  of  forces  between  four  equal  fractions  of 
the  maximum  computed  contact  force.  For  example,  if  the  maximum  contact  force  is  F  KN,  the 
narrowest  rectangle  stands  for  the  range  of  forces  between  0  and  F/4  KN,  the  next  wider 
rectangle  for  the  range  of  forces  between  F/4  and  F/2  etc.  This  notation  will  be  used  throughout 
this  report  whenever  results  from  the  discrete  element  method  are  presented. 

Figure  6  depicts  the  geometric  configuration  of  the  3-D  395-particle  array  under  o8  =  131  KPa. 
In  this  figure,  the  difference  in  color  particles  identifies  their  sizes.  The  395  sphere  3-D  array 
with  the  particles  assigned  the  properties  of  the  glass  beads  used  in  the  experiments  of  Volume  n 
(Shear  Modulus  G,  =  30  GPa,  Poisson's  ratio,  v,  =  0.15  (instead  of  0.21),  friction  coefficient  = 
0.32)  was  first  subjected  to  an  isotropic  confining  pressure  a0  =  140  KPa  and  then  the  vertical 
stress,  ov  =  cu,  was  increased  while  the  mean  stress  was  kept  constant  by  unloading  the  lateral 
(horizontal)  stresses,  an  =  o„.  This  is  analogous  to  a  laboratory  triaxial  compression  test  on  a 
glass  beads  where  the  mean  stress  is  kept  constant  by  reducing  the  lateral  stress  (cell  pressure)  <?e 
=  Oj  =  o3.  The  loading  was  stopped  at  a  octahedral  shear  strain  level  =  6,-Ej  =  0.3%.  The 
macroscopic  stress-strain  curve  resulting  from  this  simulation  appears  in  Fig.  7a,  where  the 
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octahedral  shear  stress  is  plotted  versus  the  octahedral  shear  strain.  The  volumetric  strain 
appears  in  Fig.  7b,  where  it  is  plotted  versus  the  octahedral  shear  strain.  The  same  two  plots 
include  data  from  laboratory  experiments  on  glass  beads  performed  as  part  of  this  research  and 
presented  in  Volume  II  of  this  report. 

3.2  Numerical  Simulations  with  Unloading  Investigating  Initial  and  Subsequent  Yield  Loci 

The  issue  investigated  is  the  micromechanical  mechanisms  controlling  the  movement  and 
distortion  of  the  yield  surface  of  a  granular  medium  during  loading.  Analytical  considerations 
and  laboratory  results  presented  in  Volume  II  of  this  report  indicate  that  the  yield  surface  of  a 
granular  medium  distorts  in  the  direction  of  loading  in  a  manner  analogous  to  that  of 
polycrystalline  aggregates  (metals)  discussed  earlier.  In  order  to  investigate  the  nature  these 
early  macroscopic  findings  and  to  gain  insight  into  the  micromechanical  processes  occurring  at 
the  contact  level  that  shape  the  macroscopic  behavior,  a  series  of  2-D  and  3-D  numerical 
simulations  were  performed  using  program  CONBAL  and  CONBAL-P,  on  the  477  and  531 
particle  2-D  arrays  without  rotation  and  on  a  395  particle  3-D  array  with  rotation.  The  stress 
paths  along  which  these  numerical  simulations  were  performed  are  very  similar  to  those 
proposed  by  Phillips  and  his  co-workers  and  are  shown  in  Fig.  8.  Since  the  numerical 
simulations  in  contrast  to  the  laboratory  experiments  are  strain-controlled,  all  points  were  first 
defined  in  stress  space  and  then  their  positions  were  translated  to  stress  space. 

The  "yield  criterion"  used  was  the  stress  corresponding  to  a  change  in  the  value  of  octahedral 
shear  straiv,  yM ,  of  0.024%  measured  between  the  point  from  which  the  probes  start  and  any 
point  on  the  yield  locus.  The  procedure  has  been  identical  to  that  used  in  the  laboratory 
experiments  and  described  in  detail  in  Volume  n.  The  value  of  y^  used  is  not  very  different 
from  the  value  for  the  threshold  strain,  y,,  observed  in  sand,  0.01  <yt<  0.02%,  discussed  earlier 
in  this  volume.  The  numerical  media  were  prestrained  to  a  point  in  strain  space  (point  1  in  Fig. 
8)  and  then  they  were  unloaded  to  point  2  which  was  at  a  distance  from  point  1  that  was  equal  to 
twice  the  criterion  of  octahedral  strain  used  (y*,  =  0.024%).  It  was  then  reloaded  by  0.024% 
from  point  2  to  point  O'.  Point  O'  was  used  as  the  reference  point  for  all  subsequent  small  strain 
loading  (probes)  to  establish  the  shape  of  the  subsequent  yield  locus.  The  rest  of  the  points  were 
obtained  by  starting  from  O'  and  by  probing  in  all  directions;  this  is  illustrated  in  Fig.  8  where 
the  numbers  of  the  points  are  in  the  order  used  to  obtain  those  points.  At  this  small  strain  level 
of  probing,  few  inelastic  deformations  which  are  the  result  of  geometric  changes  in  the  array  are 
expected  in  the  medium  after  probing  from  point  O'.  However,  most,  if  not  all,  nonlinearity  in 
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the  macroscopic  stress-strain  behavior  at  this  small  strain  level  is  expected  to  be  the  result  of  the 
nonlinear  force-deformation  behavior  at  the  interparticle  contacts. 

The  531  array  was  consolidated  to  132  KPa,  and  then  two  series  of  calculations  simulating  four 
tests  on  four  initially  identical  numerical  "specimens"  were  conducted  in  which  the  array  was 
prestrained  in  i)  compression  to  two  different  strain  levels,  and  ii)  in  shear  to  two  different  strain 
levels.  All  four  "specimens"  were  prestrained  first  to  7^=  0.98%  and  then  to  foa  =  2.33  %. 
The  mean  stress  was  kept  constant  to  crm=  132  KPa  throughout  each  simulation. 

Figure  9  portrays  the  results  of  the  2-D  simulations  on  the  531  particle  array  in  the  x12,  (au-a22) 
stress  space.  There  are  four  yield  surfaces  in  this  figure,  the  circular  being  the  initial  surface 
defined  by  four  separate  monotonic  simulations.  The  "subsequent"  yield  surfaces  correspond  to 
prestraining  to  two  strain  levels,  foa  =  0.98%  and  2.33%  as  marked  in  two  directions 
(compression,  on-<J22 ,  and  shear,  x12 ).  It  can  be  seen  that  there  is  a  significant  change  in  the 
shape  of  the  prestrained  yield  surface:  it  has  shortened  in  the  direction  of  loading  while  its 
dimension  in  the  other  direction  has  not  changed.  Moreover,  this  shrinking  of  the  yield  surface 
increases  with  increasing  prestraining.  Finally,  the  surface  became  flatter  in  the  direction 
opposite  to  the  direction  of  loading.  While  the  initial  surface  has  a  diameter  of  10.6  KPa,  the 
width  of  the  first  subsequent  locus  in  the  compression  direction  has  decreased  to  5.4  KPa  and 
that  of  the  second  to  4.8  KPa. 

The  477-particle  array  was  consolidated  to  om  =  91  KPa  and  then  the  same  procedure  as  in  the 
531  particle  array  was  followed.  The  array  was  prestrained  to  0.99%  then  to  =2.0%. 
The  mean  stress  was  kept  constant  to  om  =  91  KPa  throughout  the  simulation.  The  yield  surfaces 
obtained  on  the  477-equal  sphere  array  are  shown  in  Fig.  10.  Again  there  are  five  yield  surfaces 
in  this  figure,  the  circular  being  the  initial  surface  defined  by  the  four  separate  monotonic 
simulations.  Two  of  the  other  yield  surfaces  correspond  to  the  yield  loci  obtained  after  the  array 
was  prestrained  in  compression  at  ■poet  =  0.99%  and  2.0%  respectively,  while  the  other  two  to 
those  obtained  after  the  array  was  prestrained  in  shear  at  "foa  =0.98%  and  2.3%.  All  subsequent 
yield  surfaces  have  shortened  in  the  direction  of  prestraining,  with  the  width  of  the  surface 
decreasing  as  the  straining  increases.  The  surfaces  have  become  pointed  in  the  direction  of 
loading  and  flatter  in  the  opposite  direction.  The  width  of  the  surface  in  the  direction 
perpendicular  to  the  prestraining  direction  has  not  changed.  These  observations  are  in  good 
general  agreement  with  the  results  on  the  531  particle  array  (Fig.  9),  but  the  results  are  not  as 
smooth  or  symmetric.  This  should  be  attributed  to  the  fact  that  the  477  equal  particle  array  is 
not  statistically  uniform  due  to  crystallization  (Fig.  4).  This  crystallization  has  lumped  the  array 
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into  N  constituents,  where  N  is  the  number  of  different  packings  within  the  aggregate,  and  thus 
has  increased  the  characteristic  size  of  the  smallest  constituent,  from  the  diameter  of  the  particle 
to  the  dimension  of  the  crystalline  region.  As  a  result,  the  array  is  not  completely  isotropic 
under  isotropic  stress  and  a  much  larger  number  of  spheres  than  477  is  needed  in  order  to  have  a 
uniform  spatial  distribution  of  the  crystallized  regions  and  achieve  a  statistically  isotropic 
medium. 

The  395  3-D  array  was  prestrained  to  =  0.25%  and  the  subsequent  yield  locus  (surface)  was 
obtained  by  following  the  same  stress  paths  as  in  Fig.  8.  In  this  simulation  the  particles  were 
allowed  to  rotate  and  a  simplified  contact  model  was  used  to  model  the  interparticle  force 
deformation  behavior.  The  yield  surface  obtained  appears  in  Fig.  11  together  with  the  initial 
locus  which  was  obtained  from  radial  monotonic  proportional  tests.  It  can  be  seen  that  while  the 
locus  has  distorted  and  developed  an  apex  in  the  direction  of  loading  it  has  become  flatter  in  the 
opposite  direction.  Its  width  in  the  direction  perpendicular  to  prestraining  has  remained 
approximately  the  same,  indicating  that  there  is  a  minimal  cross-effect 

3.3  Micromechanical  Observations 

The  results  of  these  three  simulations  are  in  good  qualitative  agreement  with  the  laboratory  data 
presented  in  Volume  n  and  are  very  similar  to  the  yield  surfaces  obtained  in  aluminum  by 
Phillips  (Fig.  1).  As  can  be  seen  in  Fig.  1,  the  yield  locus  (surface)  of  aluminum  distorts  in  the 
direction  of  loading  by  forming  an  apex  at  the  loading  point  and  by  becoming  flatter  in  the 
opposite  direction.  The  width  of  the  surface  remains  the  same.  While  the  exact  mechanism  of 
the  yield  locus  distortion  in  metals  is  not  precisely  known,  it  is  quite  possible  that  after 
prestraining,  many  of  the  dislocations  are  pinned  at  obstacles  and  are  no  longer  mobile  in  the 
forward  strain  direction  (Helling  et  al.  1986).  The  resulting  high  strain  gradient  elevates  the 
yield  point  in  that  direction.  When  the  loading  is  reversed,  these  dislocations  become  mobile 
again,  the  yield  point  is  reached  easily  and  the  yield  locus  flattens  in  that  direction.  If  the 
loading  is  continued  along  a  path  perpendicular  to  the  prestrain  direction,  the  mobile  dislocation 
density  is  low  and  the  yield  locus  cannot  be  reached  easily.  Since  the  prestraining  has  not 
mobilized  dislocations  in  that  direction,  the  yield  point  should  be  at  approximately  the  same 
point  as  the  yield  point  on  the  stress  free  material,  and  in  this  direction  the  yield  locus  should 
remain  at  the  same  point 

The  case  of  granular  medium  is  in  a  way  similar,  although  large  strains  in  soils  are  an  order  of 
magnitude  smaller  than  in  metals.  Despite  the  fact  that  granular  media  properties  are  mean 
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stress-dependent,  this  effect  can  be  discounted  if  the  experiments/simulations  are  performed  at 
constant  mean  stress.  At  each  of  the  points  in  these  series  of  simulations  a  wealth  of  statistical 
information  has  been  accumulated  at  every  point  on  the  surface.  Figures  12,  13,  14  and  15 
present  selected  statistical  information  at  the  initial  state  before  prestraining  and  at  points  1,  2, 
and  3  respectively  on  the  2-D  yield  locus  of  Fig.  9.  This  information  includes  the  distributions 
of  i)  contact  angle,  ii)  normalized  mobilized  angle,  iii)  coordination  number  (number  of  contacts 
per  particle)  and  iv)  the  distribution  of  the  contact  force.  The  normalized  mobilized  angle  is 
the  angle  between  the  contact  force  and  contact  normal  normalized  by  the  angle  of  friction. 
Therefore,  when  the  normalized  mobilized  angle  is  equal  to  1,  particles  are  sliding,  and  if  it  is 
equal  to  0  contacts  are  subjected  only  to  a  normal  force.  The  difference  in  signs  indicates  a 
different  force  direction. 

As  can  be  seen  by  the  distribution  of  contact  angle  in  Figs.  12,  13  and  14  there  is  a  clear 
geometric  texture  (fabric  anisotropy)  developing  in  the  direction  of  prestrain  compared  to  the 
initial  state  of  the  material  (Fig.  12).  The  fabric  at  this  point  is  controlling  the  subsequent 
"yield"  behavior  of  the  material,  since  the  strain  probes  from  point  O'  defining  the  yield  locus 
(surface)  use  a  strain  low  enough  (y^  =  0.02%)  that  cannot  cause  significant  geometric  changes 
to  the  fabric  of  points  on  the  locus.  Consequently,  the  distribution  of  the  contact  angle  is  very 
similar  for  all  points  on  the  yield  surface  (Figs.  13, 14, 15). 

The  distribution  of  the  normalized  mobilized  angle  is  very  different  from  point  to  point.  While 
the  number  of  spheres  sliding  does  not  change  significantly  (both  the  coordination  number  and 
the  number  of  normalized  mobilized  planes  equal  to  1.0  remains  the  same),  what  changes 
dramatically  is  the  distribution  of  both  the  mobilized  angle  and  contact  force  that  controls  the 
contacts  that  will  be  sliding.  For  example,  when  the  sample  was  loaded  in  compression,  contacts 
were  created  in  the  direction  of  the  major  applied  principal  stress  while  they  were  lost  in  the 
other  directions.  The  Distribution  of  the  mobilized  angle  and  contact  force  changed  accordingly. 
Consequently,  when  point  1  was  reached  and  an  anisotropic  fabric  has  been  created.  Once  the 
loading  is  reversed  and  the  force  redistribution  takes  place,  it  is  easier  for  the  sample  to  yield  in 
the  direction  opposite  to  the  loading,  due  to  the  texturing  which  has  occurred  in  the  direction  of 
loading.  This  can  be  observed  through  the  geometry  of  the  numerical  sample  at  both  points 
(Figs  16  and  17). 

The  distributions  of  the  normalized  mobilized  angle  at  points  1  and  2  are  very  similar  to  being 
the  mirror  image  of  the  other.  This  suggests  that  the  forces  at  the  particle  contacts  have  changed 
sign  or  reversed  direction.  Moreover,  the  fact  that  the  specimen  geometry  has  not  changed 
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significantly  between  these  two  points,  indicates  that  there  are  no  (or  very  few)  spheres  that  slid 
between  points  1  and  2.  Therefore,  there  is  slipping  (or  rubbing),  but  not  sliding,  at  the 
interparticle  contacts  along  the  same  "paths"  that  sliding  occurred  during  loading  to  point  l. 
Since  there  are  no  pinned  sliding  systems,  there  are  no  deformation  gradients  and  the  yield  locus 
flattens. 

When  the  sample  is  loaded  by  a  small  strain  shear  probe,  the  mobilized  angle  distribution  is 
analogous  to  that  of  the  sample  before  prestraining.  This,  coupled  with  the  fact  that  this 
direction  of  loading  is  orthogonal  to  the  prestraining  direction  no  sliding  systems  have  been 
activated  in  this  direction,  results  in  a  yield  point  similar  to  that  of  the  original  sample. 
Consequently,  the  size  of  the  surface  does  not  change  perpendicular  to  the  direction  of 
prestraining. 

The  same  is  true  for  the  477  equal  particle  array  as  can  be  seen  from  the  analysis  of  the 
geometrical  and  statistical  information  presented  in  Figs.  16-20  for  points  1,  2,  and  3  on  the 
yield  locus  of  Fig.  9. 

The  nonlinear  distinct  element  method  has  been  very  useful  in  predicting,  as  well  as  interpreting 
aspects  of  the  behavior  of  granular  medium  which  could  not  be  determined  otherwise.  They 
could  not  be  determined  by  a  distinct  element  procedure  using  a  linear,  non-pressure  dependent 
force  deformation  law  at  the  interparticle  contacts,  since  this  program  would  be  unable  to 
capture  the  redistribution  of  contact  forces  which  occur  during  loading.  Nor  could  these 
simulations  be  performed  without  a  supercomputer,  since  the  calculation-demanding  relaxation 
scheme  used  in  the  distinct  element  method,  burdened  by  a  nonlinear  plasticity  program  at  every 
contact,  could  not  be  solved  on  a  regular  computer.  As  an  indication  of  the  complexity  of  the 
computations,  the  CPU  time  needed  to  construct  Fig.  8  or  Fig.  9  was  approximately  30 
CPU-hours  on  the  IBM  3090-600E.  Therefore,  as  previously  hypothesized  by  the  authors,  the 
nonlinear  behavior  of  the  granular  medium  is  a  result  of  the  particulate  nature  of  the  medium  and 
can  not  be  interpreted  or  reproduced  analytically  unless  this  particulate  nature  is  fully  modelled 
and  taken  into  account 
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4.  CONCLUSION 

Nonlinear  discrete  element  simulations  have  been  used  to  provide  an  insight  on  the  nonlinear 
modelling  of  granular  soil.  These  simulations  were  based  on  an  incremental  solution  to  the 
nonlinear  problem  of  two  spheres  in  contact  (program  CONTACT),  incorporated  into  discrete 
element  TRUBAL  which  was  further  optimized  for  vector  and  parallel  processing  on  the  IBM 
3090  supercomputer.  It  has  been  found  that  this  approach  not  only  interprets  successfully  the 
nonlinear  behavior  of  soil,  but  also  provides  a  wealth  of  information  on  the  fabric  changes 
during  loading.  The  yield  surface  of  a  granular  medium,  needed  for  defining  the  constitutive 
relation  of  such  a  medium,  distorts  by  forming  an  apex  in  the  direction  of  loading  while 
becoming  flatter  in  the  opposite  direction.  This  is  contrary  to  the  practice  followed  in  modelling 
granular  media  where  the  yield  surface  of  soils  are  typically  assumed  to  retain  the  same  shape. 
The  origin  of  this  distortion  phenomenon  lies  in  the  texturing  (or  fabric  anisotropy)  which  occurs 
in  the  direction  of  prestraining,  as  well  as  in  the  redistribution  of  interparticle  contact  forces  in 
the  absence  of  significant  particle  movement  during  the  small  strain  probes  needed  to  define  the 
yield  surface.  These  phenomena  cause  certain  slip  systems  to  be  activated  which  in  turn  produce 
the  characteristic  apex  which  appears  in  the  yield  surface  in  the  loading  direction.  Therefore,  the 
distribution  and  magnitude  of  the  contact  forces  are  critical  for  a  good  understanding  of  the 
macroscopic  response  of  the  medium.  Accurate  modelling  of  the  contact  force  distribution  can 
be  achieved  only  if  the  behavior  at  the  contact  is  fully  understood  and  rigorously  modelled. 
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Fig.  la.  Loading  paths  for  (a)  initial,  and  (b)  subsequent  yield  surfaces  (Phillips  1968) 


Fig.  lb.  Experimentally  obtained  initial  and  subsequent  yield  surfaces  for  aluminum  at  different 
temperatures. 
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particle 

Acceleration  and  Spin  Computed  as  if  Neighboring 
Particles  did  not  exist.  Time  step  is  assumed  to  be 
small  such  that  accelerations  and  velocities  are 
constant  in  this  small  time. 

fiffif  t  t 

•  Particle  occupies  new  position  due  to  acceleration 

•  New  Pi's,  Mi’s  are  computed  for  new  positions 

•  IP1, 1  Ml  and  process  is  repeated. 


Figure  3. 


Main  Features  of  the  Discrete  Element  Method 
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Fig.  4.  2-D  random  array  of  477  equal,  elastic,  rough  spheres  assigned  the  properties  of  quartz, 
subjected  to  isotropic  compression.  This  figure  represents  the  "element"  of  the  periodic  space. 
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Fig.  5.  2-D  random  array  of  531  elastic,  rough  spheres  with  two  radii  (Rt  /  R*  =  1.5)  assigned 
the  properties  of  quartz,  subjected  to  isotropic  compression.  This  figure  represents  the  "element" 
of  the  periodic  space. 


Fig.  6.  Initial  configuration  of  the  3-D  random  array  of  365  elastic. 
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Octahedral  Shear  Strain  (%) 


Pig.  7a.  Comparison  between  the  octahedral  stress  - 
strain  enrvas  obtained  by  the  3-0  numerical 
simulation  and  laboratory  experiments. 


Pig.  7b.  Comparison  between  the  volumetric 

•train  *  octahedral  shear  strain  curves 
obtained  by  the  3-D  numerical 
■Emulation  and  laboratory  experiments. 
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Fig.  8.  Stress  paths  in  stress  space  used  in  determining  the  initial  and  subsequent  yield  loci  in 
the  numerical  simulations. 
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Fig.  9.  Initial  and  subsequent  yield  loci  obtained  for  the  531  particle  2-D  array.  Prestraining  as 
indicated. 
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Fig.  11.  Initial  and  subsequent  yield  loci  obtained  for  the  36S  particle  3-D  array.  Prestraining 
indicated. 
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Fig.  12.  Statistical  Information  regarding  the  state  of  the  531  particle  array  of  Fig.  5  under 
isotropic  pressure,  on  =  132  KPa. 
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Fig.  14.  Statistical  Information  regarding  the  state  of  the  531  particle  array  of  Fig.  5  at  point  2 
on  the  yield  surface  of  Fig.  9. 
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Fig.  15.  Statistical  Information  regarding  the  state  of  the  531  particle  array  of  Fig.  5  at  point  3 
on  the  yield  surface  of  Fig.  9. 
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Fig.  18.  Statistical  Information  regarding  the  state  of  the  477  particle  array  of  Fig.  4  at  point  1 
on  the  yield  surface  of  Fig.  10. 
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Fig.  19.  Statistical  Information  regarding  the  state  of  the  477  particle  array  of  Fig.  4  at  point  2 
on  the  yield  surface  of  Fig.  10. 
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Fig.  20.  Statistical  Information  regarding  the  state  of  the  477  particle  array  of  Fig.  4  at  point  3 
on  the  yield  surface  of  Fig.  10. 
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Fig.  21.  Geometry  and  contact  force  distribution  of  the  477  particle  array  of  Fig.  5  at  point  1  on 
the  yield  surface  of  Fig.  10. 


Fig.  22.  Geometry  and  contact  force  distribution  of  the  477  particle  array  of  Fig.  5  at  point  2  on 
the  yield  surface  of  Fig.  10. 


